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ABSTRACT 

We investigate the structure and the stabilities of a protoplanetary disk, which is 
heated by viscous process in itself and by its central star. The disk is set to rotate 
with the Keplerian velocity and has the surface density distribution of the minimum 
mass solar nebula. We assume the vertical hydrostatic equilibrium and the radiative 
equilibrium at each point, and solve the two-dimensional radiative transfer equation 
by means of the Short Characteristic method in the spherical coordinate in order to 
determine the disk structure. Our calculation shows that at the outer region of the disk 
with a distance from the central star of x > lAU the radiative heating from the inner 
disk dominates the viscous heating even near the midplane. It is because of the high 
temperature distribution in the optically thin surface layer and the relatively high disk 
height {zoo ~ 0.7x at x ~ lAU) as a consequence of the irradiation from the inner hot 
region of the disk. In addition, we examine the convective and the magnetorotational 
instabilities of the disk. As a result, the whole disk is convectively stable since the 
dusty region is not heated by the viscous dissipation from the midplane but by the 
radial radiative heating. On the other hand, almost all the disk is magnetorotationally 
unstable except for the region near the equatorial plane of 2AU < x < lOAU. Finally we 
discuss the growth and the size distribution of dust particles in the disk, which suggests 
that there exist cm-sized particles in the surface layer, namely, in the exposed region of 
the disk. 

Subject headings: accretion, accretion disks — circumstellar matter — instabilities — 
planetary systems: protoplanetary disks — radiative transfer 


1. INTRODUCTION 

It has been established that there exist circumstellar disks around a large fraction of young 
stellar objects (YSOs), some of which will evolve to form planetary systems analogous to our Solar 
System (Wilking, Lada, & Young 1989; Stauffer et al. 1994; Haisch, Lada &: Lada 2001). Infrared 
excesses over the photospheric emissions, often observed in the spectral energy distributions (SEDs) 
of YSOs, are considered to be radiated from dusty circumstellar disks (Cohen 1974; Adams & Shu 
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1985; Kenyon & Hartmann 1987; Beckwith & Sargent 1993; Hartmann 1998). Furthermore, the 
disk-structure images with radii of ~ lOOAU are sometimes observed around YSOs at optical and 
near-infrared wavelengths as silhouettes against the background bright nebula (O’Dell et al. 1993; 
O’Dell 1998; Bally, O’Dell & McCaurean 2000) or as optically thick dust lanes with optically thin 
reflection nebulae, absorbing and scattering light from their central stars (Burrows et al. 1996; 
Padgett et al. 1999). 

These circumstellar disks play vital roles on the star and planetary system formation: for 
example, they transform angular momentum outward and accrete mass to their central stars (e.g., 
Spitzer 1978; Shu, Adams, & Lizano 1987; Bodenheimer 1995; Papaloizou & Lin 1995; Nomura & 
Mineshige 2000; Stone et al. 2000), and they form planets in them through the various processes 
of dust grain growth (e.g., Safronov 1969; Weidenschilling & Cuzzi 1993; Beckwith, Henning & 
Nakagawa 2000). The turbulence in the disks is one of the keys to elucidate their structures 
and evolutions. The turbulent viscosity in differentially rotating disks transforms the angular 
momentum and mass, and releases the gravitational energy of accreted mass as heat through the 
viscous dissipation (e.g., Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; Pringle 1981). A 
number of mechanisms to excite turbulent motion have been investigated for various regions of disks 
in varied evolutional stages, such as the thermal convection, the magnetorotational instability, the 
shear motion due to differential rotation, the infalling envelope, the stellar wind, and so forth (e.g., 
Lin & Papaloizou 1980; Cabot et al. 1987; Balbus & Hawley 1991; Sano et al. 2000; Zel’dovich 
1981; Sekiya 2000; Elmegreen 1978; Cameron 1978), but we do not yet have a universally accepted 
scenario. Meanwhile, many studies have been done on the structure of such optically thick viscous 
(active) protostellar disk under the assumptions of the vertical hydrostatic equilibrium and the 
energy balance between the viscous heating and the radiative cooling, sometimes in the context 
of explaining the FU Ori outburst phenomena (e.g., Lin &: Papaloizou 1980; Pringle 1981; Lin & 
Papaloizou 1985; Kawazoe & Mineshige 1993; Bell & Lin 1993). 

In order to reproduce the structure of protoplanetary disk, however, it is not sufficient to 
evaluate only the vertical viscous heating from the midplane since the disks are irradiated from 
the central star and the inner hot region of the disk (e.g., Kenyon & Hartmann 1987; Chiang & 
Goldreich 1997; D’Alessio et al. 1998; Bell 1999). This irradiation process is noteworthy since it 
heats up and flares up the disk surface, reproducing the observed flat SEDs of far-infrared radiation 
from YSOs (the radiation flux of (xv^), which is not recreated by radiation from geometrically 
thin accretion disks (e.g., Lynden-Bell & Pringle 1974; Adams, Lada, & Shu 1987). The flared 
features of optically thick disks are in fact observed as shadows against scattering starlight in 
optically thin circumstellar matter (e.g., Burrows et al. 1996; Padgett et al. 1999). The profiles 
of such passive disks irradiated only from their central stars have been investigated under the 
assumption of the vertical hydrostatic equilibrium and with a simple model of radiative transfer of 
starlight, without and with the turbulent viscous heating (e.g., Kenyon & Hartmann 1987; Chiang 
& Goldreich 1997; Hubeny 1991; D’Alessio et al. 1998). 

In this paper we strictly treat the radiation processes by solving the two-dimensional radiative 
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transfer equation so as to investigate the structure and the instabilities of an axisymmetric disk 
under the influences of the viscous heating and the reprocessing of radiation from the central star 
and the inner hot region of the disk. In the next section we present the disk model, which is heated 
by the viscous dissipation and its central star, and assumed to satisfy the vertical hydrostatic 
equilibrium and the local radiative equilibrium. In §3 we calculate the density and temperature 
distributions of the disk. The convective and the magnetorotational instabilities of the disk, which 
will excite turbulent motion, are investigated in §4, and the dust size distribution in the disk 
is discussed by means of the turbulent-eddy-trapping model for dust growth in §5. Finally we 
summarize the results in S6. 


2. DISK MODEL 

We consider an axisymmetric disk surrounding a central star with physical parameters of the 
typical T Tauri stars; the mass of M* = O.SMq, the radius of i?* = 2Rq, and the temperature 
of = 4000K (e.g., Kenyon & Hartmann 1995). We investigate the density and temperature 
distributions of the disk under the assumptions of the hydrostatic equilibrium in the vertical di¬ 
rection (§2.1) and the radiative equilibrium at each point (§2.3). We also hypothesize that dust 
and gas have the same temperature distribution allover the disk in this paper. For a matter of 
convenience, we adopt the cylindrical coordinate (x,z) in §2.1 and the spherical coordinate (r, 0) 
[where r = (x^ -|- and 0 = {x/z)] in §2.3 , in both of which we put the central star at 

the origin. 


2.1. Hydrostatic Equilibrium in the Vertical Direction 


In order to determine the density distribution, p{x,z), we assume the vertical hydrostatic 
equilibrium. 


dP dp 
dz p dz 


( 1 ) 


The sound speed Cg is defined as Cg^ = dP/dp = kT/m^, where P, T, k and represent the 
pressure, the temperature, the Boltzmann’s constant, and the mean molecular mass, respectively, 
and we adopt = 2.3mH (mn is the hydrogen mass). The circumstellar disks around T Tauri 
stars are considered to be rotationally supported. Then the gravitational force in the vertical 
direction is written as Pz = with the Kepler frequency. 


H(x) = (GM*/x3)^/2 _ 14 X io-7s-1(x/1AU)-3/2^ 


( 2 ) 


where we adopt M* = 0.5Mq for the mass of the central star. For the boundary condition, we 
put p(x, Zoo) = 3.8 X 10“^^ g cm“^, which corresponds to the typical number density of molecular 
clouds, n ~ 10^ cm“^ (where the mean molecular mass of = 2.3mH is assumed). The height of 
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the disk Zoo is determined by the condition, 

/ ^oo 

p(x, z)dz = E(x), (3) 

'^oo 

where the surface density of the disk S is put to agree with that of the minimum mass solar nebula 
(e.g., Hayashi, Nakazawa, & Nakagawa 1985), 

S(x) = 1.7 X 10^ g cm-2(x/lAU)-3/2. (4) 

Solving equation (1) with the above conditions and the temperature distribution T{x,z), which is 
derived in subsection §2.3, by means of the fourth-order Runge-Kutta method, we get the density 
distribution p{x,z). 


2.2. Heating Sources 

For the heating sources of the disk we consider two kinds of processes; the viscous dissipation in 
itself and the gravitational energy release associated with the contraction of its central T Tauri star. 
In protoplanetary disks it is believed that some kinds of instabilities (see §4) cause the turbulent 
motion and then induce the angular momentum and mass transfer. This process is accompanied 
with the heating of the disks via viscous dissipation of turbulence. Following the so-called a-viscous 
model, we represent it by putting heating source at the midplane of the disk (at z = 0,0 = 7r/2) 
with a heating rate of 

Q+, = (9/4)SaCsgfl, (5) 

where Cso is the sound speed at the midplane. In this model the kinetic viscosity is prescribed as 
zzvis = 0 (CsqH = acso/f^, where H is the scale height of the disk (e.g., Shakura & Sunyaev 1973; 
Pringle 1981; Kato, Fukue & Mineshige 1998), and a = 10“^ is adopted here. The vicinity of the 
midplane at the inner region of the disk is mainly heated by diffusive radiation from this viscous 
heating source as we will see in §3. Besides, we put thermal radiation source at the stellar surface 
(at r = ii*) with a heating rate of 

Qtt.r = <yTt, ( 6 ) 

where a is the Stehne-Boltzmann constant, and examine the reprocessing of radiation from the 
central star in the disk by solving the two-dimensional radiative transfer equation as we describe in 
the following subsection (§2.3) and Appendix. For the stellar temperature T* = 4000K is adopted. 
We here comment that the irradiation not only from the central star but also from the inner hot 
disk on the surface region of the outer disk is automatically simulated in our calculation, and the 
latter dominates the former in our model as we will see in §3. 
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2.3. Radiative Equilibrium 


Evaluating the temperature distribution T(r, 0), we assume the radiative equilibrium (that 
is, the emitted and absorbed radiation are balanced) at each point (r, 0) [r = is the 

distance from a central star and 0 = tan“^(x/ 2 ;) is the angle from the z-axis] as 

POD POC 

dvr/ K,u{r,@)Bt,[T{r,Q)]di/ = / K,u{r,Q)Ju{r,Q)di/, (7) 

Jo Jo 

where i' is the photon frequency, By[T) the Plank function for blackbody radiation, and Ky the 
monochromatic opacity given by the model in the following subsection (§2.4). In equation (7) local 
thermodynamics equilibrium riy{r, 0) = Ku{r, Q)By[T{r, 0)] is assumed, where r]y is the monochro¬ 
matic emissivity (e.g., Mihalas 1978). The mean intensity Jy{r,Q) is given by integrating the 
specific intensity Iy{r, 0; ;U, 4>) of radiation field for all solid angles diu = dfidcj) as 

p2'k p\ 

Jy{r,e) = —J^ J ^Iy{r,e]n,J))dfid(j), (8) 

where // is the cosine of the angle between the ray direction and the vector Op {O = [0,0],P = 
[r, 0]), and (j) is the angle between the projection of the ray to the plane whose normal is OP and 
the line ((/> = 0) on the plane parallel to the midplane of the disk (z = 0, 0 = 7r/2) (cf. Fig. 3 of 
Dullemond & Troulla 2000). We here divide the total specific intensity ly into three parts as 

= (9) 


where I* and I™ represent radiation directly from the heating sources described in §2.2, the central 
star and the viscous dissipation at the midplane of the disk, respectively, and arises from thermal 
emission of dust grains in the disk (e.g., Shu 1991; Efstathiou & Rowan-Robinson 1990). We note 
that we neglect the radiative scattering for the simplicity in this paper. Each of them is derived 
from radiative transfer equation. 


d^ 

ds 


pKy{Sl 


It) 


( 10 ) 


where the source functions S* = S'™ = 0 and S*'^ = By. From the above equation, the intensity I* 
and ly^^ are derived with = T{r = R^,,Q) and Tvis(x) = T(r, 0 = 7r/2) as 






( 11 ) 


for jx < —(r^ — [that is, (/U,(/>) originates from the central star] and Iy{r,@; p, 4>) = 0 

otherwise, and 

/r(r,0;//,</.) = R,[r™(r)]e-"‘'(^’®) (12) 

when {iJ,,4>) originates form the midplane (r, 0 = tt/2) = {x,z = 0) and /™(r, 0;//, i;A) = 0 other¬ 
wise. The viscous heating temperature rvis(?’) is related with Cso in equation (5) as Csg = kT^is/m^. 
The intensity is derived by integrating equation (10) along a ray, s (which is directed to [p, (/>]), 
as 

0; p, </>) = r Ky{r\ Q')p{r', Q')By[T{r', ds'. (13) 

Jo 
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In equations (11)-(13) Tu(r,0) represents the specific optical depth at the point (r, 0) from the 
radiation source (r',0'), given by 

T^(r,e)= [ K^{r',Q')p{r',Q')ds\ (14) 

Jo 

Together with the density distribution p{r,Q) obtained in §2.1, we solve equations (7)-(14) self- 
consistently to get the temperature distribution T{r, 0) in the disk. The numerical method to solve 
equation (7) is described in Appendix A.l, and that to calculate equations (11)-(14) is in A.2. We 
note that we do not consider energy transport by convection here, but it turns out in §4.1 that this 
is permissible since the whole of the disk is in fact convectively stable in our disk model. 


2.4. Opacity 

In this paper we consider that the disk consists of dusty medium (in which dust is well 
mixed with gas) at low temperature (T < 2300K) and only of gaseous matter at T > 2300K. 
For dusty medium we use the frequency dependent absorption coefficient Ku of Adams & Shu 
(1985, 1986; see also Draine & Lee 1984). In their model the dust is constituted by three 
types of components of graphite and olivine silicate grains, and water icy mantle which coats 
the grains. They are assumed to evaporate at the characteristic dust destruction temperatures 
of 2300 K (graphite), 1500 K (silicate), and 150 K (ice). For gaseous matter we provision¬ 
ally take the weighted transmission averaged opacity, the Rosseland mean opacity kr, defined as 
1 /kr = J^{l/Ky){dBy/dT)dv/ {dBy/dT)dv. The density and temperature dependent model 
of Bell & Lin (1994), kr = Kip°'T^, is adopted here, where k*, a and b are given for each of eight 
states of gas. It represents the Alexander/Cox/Stewart opacity and covers the wide temperature 
range of lO^K < T < lO^K (see Bell &: Lin 1994 in details and see also Lin &: Papaloizou 1985). 
The applicability of this approximate treatment of the opacity for gas is discussed in §3. 


3. DENSITY AND TEMPERATURE PROFILES OF THE DISK 

Solving the equations of the vertical hydrostatic equilibrium (1) and the local radiative equi¬ 
librium (7) iteratively together with the two-dimensional radiative transfer equation (10) as we 
mentioned in the previous section and Appendix, we obtain the self-consistent density and temper¬ 
ature distributions of the disk heated by the viscous process in itself and by its central star. Figure 
la shows the contour plots of the resultant temperature {solid lines) and density {dotted lines) pro¬ 
files in the x-z plane of the ‘irradiated disk’ [which means that it is the result of the two-dimensional 
calculation of radiative transfer]. The contour levels are T = 10,10^,10^,10^, and lO^K for the 
temperature, and p = 10“^®, 10“^®, 10“^^, 10“^^, 10“^*^, and 10“®g cm“^ for the density. We also 
plot the vertical temperature and density distributions at x = 0.046,0.1, 0.21, 0.46,1.0, 2.1,4.6,10, 
and 21AU in Figure 2a and 3a, respectively. We here mention the applicability of our opacity 
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model for gas (T > 2300K). Whereas we use the Rosseland mean opacity in this paper as we stated 
in §2.4, the Planck mean opacity is to be used in equation (13) and the left side of equation (7), 
and the energy-weighted one in the right side of equation (7) (see Hubeny 1990). This discrepant 
treatment of the opacity will modify the temperature profile at the optically thin region of the 
inner disk because the difference among these opacities goes up to several orders of magnitude at 
the maximum (e.g., Alexander & Ferguson 1994). The temperature near the midplane is, however, 
hardly affected since it is determined by equation (Al) as we stated in Appendix. For comparison 
we display those of the ‘non-irradiated disk’ [which is derived from the energy transport equation 
only in the vertical direction (Al)] in Figure lb, 2b, and 3b in the same manner as Figure la-3a. 
The constant temperature near the surface in Figure lb and 2b is due to our inappropriate use of 
the diffusion approximation (eq. [Al]) even at the optically thin region. 

The remarkable differences between the ‘irradiated’ and ‘non-irradiated’ cases appear in the 
rises of temperature at the surface region and at the equatorial plane outside x ~ lAU of the 
‘irradiated disk’. As a result, the vertical height of the ‘irradiated disk’ becomes higher than 
that of the ‘non-irradiated disk’. The effect of this flaring is due to the irradiation from the 
inner region and conspicuous at the outer disk, as is modeled in Kenyon & Hartmann (1987), 
especially at the optically thin surface layer (Chiang & Goldreich 1997; D’Alessio et al. 1998). 
The main source of the irradiation is, however, the inner hot disk rather than the central star 
in this case, contrasting with the previous works. We next discuss the reason for the radiative 
heating from the inner disk to dominate the viscous heating at x > lAU. Now the energy density 
due to the viscous heating at the midplane of the disk decreases with radius approximately as 
Kvis = oc which is derived from the equations of the energy balance and the radiative 

transfer in the vertical direction, {9/A)TiaCs‘^Q. = cjTJj ~ where we use S oc , 

12 oc CsQ oc Tvis, and kr oc (dust grains mainly contribute to the opacity in this 

region). Next let us consider the radial dependence of the energy density due to the radiative 
heating from the inner disk. If we postulate the extreme case that the height was large enough 
to adopt the plane-parallel approximation in the radial direction, the dependence would became 
-E'rad = (STrad/dc) f KRpdr OC r where F].ad(= const.) is the radial radiative flux and 

we use kr oc and p oc from our calculation. In this case the radiative heating would 

surely dominate the viscous heating. On the contrary, if the disk is geometrically thin enough, the 
radiative heating from the inner disk becomes negligible. Thus one of the reasons for our result 
is the high vertical height of the disk. In addition the temperature at the surface layer is higher 
than that at the equatorial plane, which allows the effect of radiative heating from the inner region 
stronger. 

In conclusion our results show that the surface layer of the disk is heated up and flared via the 
irradiation from the inner hot disk, and that the radiative heating from the inner disk dominates 
the viscous heating at x > lAU because of the high disk height and the high temperature at the 
disk surface in this case. 



4. DISK INSTABILITIES 


In protoplanetary disks it is believed that some kinds of instabilities cause the turbulent motion, 
and then, induce the angular momentum and mass transfer, and/or the energy transport. A 
number of processes have been proposed, among which we examine in this section two of intrinsical 
mechanisms in the disks before the dust settling stage: the convective instability (§4.1) and the 
magnetorotational instability (§4.2), making use of the temperature and density profiles obtained 
in the previous section. 


4.1. Convective Instability 

It has been suggested that the protoplanetary disks are convectively unstable because they 
consist of dusty components whose opacity has a temperature dependence of kr oc T^, /? > 1 , and 
then causes a superadiabatic gradient against the direction of gravitational force (Lin & Papaloizou 
1980, see also Cameron &: Pine 1973; Ruden, Papaloizou, & Lin 1988; Kley, Papaloizou, & Lin 
1993; Stone & Balbus 1996). Now we examine the convective instability of the disk, which has the 
temperature and density distributions shown in §3, using the criterion of 

A,-A,d>0, = = (15) 

a log P 7 

where 7 is the specific heat ratio and the pressure satisfies P = 'ycg^p (e.g., Cox 1980). We here 
take 7 = 7/5 since the main component of the disk is molecular hydrogen. As a result, we find 
that the whole of the ‘irradiated disk’ [obtained from the two-dimensional calculation of radiative 
transfer] is convectively stable, while the dusty region (T < 2300K) of the ‘non-irradiated disk’ 
[derived from only the vertical radiative transfer] is unstable (Fig. 4b, dotted stripe region). The 
unstable region disappears in the ‘irradiated disk’ because its dusty region is not heated from the 
midplane but from the inner region as we see in §3. 


4.2. Magnetorotational Instability 

Magnetorotational instability is one of the promising turbulent sources in accretion disks (Vel- 
hikov 1959; Balbus & Hawley 1991; Balbus & Hawley 1998 and references therein). When a 
differentially rotating, fully ionized, magnetized disk satisfies the condition, 

27ruA/H < H, (16) 

the disk is magnetorotationally unstable, where H and H are the rotational frequency and the disk 
scale height, respectively, and va is the Alfven speed, defined with the magnetic field B and the 
density p as ua = BKAirp)^^'^ (Balbus & Hawley 1991). In protoplanetary disks, on the other hand, 
the ionization degree of material is low enough (e.g., Umebayashi & Nakano 1988) that the effects 
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of the ohmic dissipation and the ambipolar diffusion should be taken into account (e.g., Blaes & 
Balbus 1994; Jin 1996; Gammie 1996; Sano et al. 2000). Here we examine the magnetorotationally 
unstable region in the disk with the temperature and density distributions obtained in §3 under 
the assumption that the ratio Csq/'^a = 10 allover the disk. According to Sano et al. (2000), in 
which the stabilizing effect of the ohmic dissipation in protoplanetary disks is investigated, the disk 
is magnetorot ationally unstable when 

27r max[uA/H, ry/uA] < H, (17) 

where r] is the magnetic diffusivity. We here consider only the collisions between electrons and 
neutrals, mainly hydrogen molecules and helium atoms, as contributor to the diffusivity r] [the 
contributions of ions and charged grains are negligible at x > lAU (Sano et al. 2000)]. Thus r/ is 
given by 

rj = c^/dTTfTc = {c^me<(JV>ej4:Tre^){njne), (18) 

where c, cTc, m-e, e, n, and Ue are the light speed, the electrical conductivity, the electron mass and 
charge, the number densities of neutrals and electrons, respectively. For the momentum transfer 
rate between electrons and neutrals, <av>e, we use the experimental formula in Appendix of Sano 
et al. (2000; see also Hayashi 1981). Meanwhile, we evaluate the ionization degree ^ = Ue/n at 
T < lO^K with the ionization rate C; the temperature T, and the molecular number density n as 
^ = 3.4 X (Gammie 1996), where C is represented as C = CcRexp(—rcR)+CR- Each 

of CcR and Cr is the ionization rate by cosmic rays and radioactive elements, respectively. Following 
Umebayashi & Nakano (1980), we adopt CcR = 10“^^s“^ and Cr = 6.9 x 10“^^s“^. The ionization 
depth by cosmic rays, tcr, is written as tqr = p(x, z)dz/xcR, where XCR is the attenuation 
length and we use xcR = 96g cm“^. At T > lO^K, where the thermal collision dominates the 
cosmic rays and the radioactive elements on ionization, we adopt the thermal ionization degree ^ of 
Umebayashi (1983). In addition, in order to see the effect of ambipolar diffusion, we examine the 
ratio of the collision frequency between ions and neutrals to the epicyclic frequency, XdPi/H, where 
7d ~ 3.5 X 10^^ cm^ g“^ s“^ is the drag coefficient (e.g., Shu 1983) and p\ is the ion density (Blaes 
& Balbus 1994; Hawley & Stone 1998). 

In Figure 4a we plot the resultant magnetorotationally unstable region (solid stripe) in the 
x-z plane of the ‘irradiated disk’, whose density and temperature distributions are obtained from 
our two-dimensional numerical calculation. The figure shows that almost all the disk is magne¬ 
torotationally unstable except for the region near the equatorial plane of 2AU < r < lOAU. The 
ionization degree is high enough for the disk to be unstable in the inner region because of the 
thermal collisions, and in the surface and outer regions owing to the cosmic rays. Equation (17) 
suggests that the stable region spreads as Csq/va becomes larger and vice versa (although we here 
examine only the case of Cso/va = 10) since the relation ua/H < p/va is satisfied allover the disk in 
this case (see also Sano et al. 2000). The stabilization due to the ambipolar diffusion is ignorable 
in this case (xdPi/H < 10“^ in allover the region which the ohmic dissipation dose not stabilize; see 
Stone & Hawley 1998). Eor comparison we also plot the unstable region of the ‘non-irradiated disk’, 
which is obtained from the energy transport equation only in the vertical direction, in Eigure 4b. 
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The stable region in the ‘irradiated disk’ is smaller than that in the ‘non-irradiated disk’ because 
of its hotter temperature profile. We note that our assumption that the disk is heated by viscous 
dissipation at every radius (§2.2) is conflict with the fact that the stable region exists in 2AU < r < 
lOAU. But our results will not be affected by this inconsistency since this region is heated not 
viscously but radiatively as we show in §3. 

In summary we conclude that the turbulent source in a protoplanetary disk is not convective 
instability but magnetorotational instability under the conditions that the disk has the surface 
density distribution of minimum mass solar nebula and the viscous parameter of a = 10“^. 


5. DUST SIZE DISTRIBUTION 

Protoplanetary disks are considered to consist at first of a few x 10“^;um sized dust grains, orig¬ 
inating from the interstellar medium, which grow larger particles and at last lead to the planet for¬ 
mation through the processes of sticking, gravitational attraction, and gas accretion (e.g., Safronov 
1969; Weidenschilling & Cuzzi 1993; Beckwith, Henning & Nakagawa 2000 and references therein). 
For the dust growth process numerous theoretical models have been proposed, for example, the 
collisional models in the laminar nebula and in the turbulent accretion disk, and the particle 
concentration model owing to turbulent motion (e.g., Kusaka, Nakano, & Hayashi 1970; Weiden¬ 
schilling 1980; Nakagawa, Nakazawa, & Hayashi 1981; Mizuno 1989; Squires & Eaton 1991; Klahr 
& Henning 1997; Supulver & Lin 2000). On the other hand, the evidences of dust grain growth in 
protoplanetary disks are observationally presented by means of the changes in their spectral energy 
distributions (e.g., Beckwith & Sargent 1991; Beckwith et al. 2000; D’Alessio, Calvet, & Hartmann 
2001; Throop et al. 2001). In this section we discuss the dust size distribution in the disk, which 
has the density and temperature structures in §3 and the turbulent region obtained in §4 for the 
‘irradiated disk’, dealing with the turbulent-eddy-trapping model of Klahr & Henning (1997). 

Dust particles can be trapped in a turbulent eddy in an accretion disk when the dust-gas 
friction time Tf satisfies the condition, 

Tf </wmin(l/(/;j, 1/2DAU), (19) 

where I and uo are the size and the rotational frequency of the turbulent eddy, respectively, and D 
is the orbital frequency around the central star, is the vertical gravitational force and 

DAK is the drag force caused by the difference in orbital velocity, AV, between gaseous matter, 
undergoing the radial gas pressure, and the Keplerian rotating dust particles (e.g., Weidenschilling 
1977). The dust-gas friction time Tf is given by 

A = apPd/csP, (20) 

where Up, pd, Cg, and p are the radius and the density of dust particle, the thermal velocity and 
the density of gas, respectively. From equations (19) and (20) we can estimate the maximum size 
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of dust grains, that can be trapped in turbulent eddies in the disk, as 

®p,max — ruin( 1 /^ 2 ,, l/2nAV^), (21) 

where we take I = z^o (disk height), oj = fi, and pd = Ig cm“^. In Figure 5 we plot the contour line 
of the maximum dust size distribution in the x-z plane of the ‘irradiated disk’, which is estimated 
from equation (21) for the turbulent (magnetorotationally unstable) and dusty (with temperature 
of T < 2300K) region of the disk. The contour levels are Op^max = 10“^, 10“^, 1,10, and lO^cm 
and the thick line displays z = Zoc- The figure shows that dust particles can grow to be cm-size in 
the surface layer, which is because both of the temperature and the disk height are high enough 
there. These particles are expected to be under the influences of the outside activities, such as 
x-ray radiations, the FU Ori outbursts, and so on, which may lead to the chondrule and/or CAI 
formations (e.g., Jones et al. 2000). The dust size growth will also be connected with the opacity, 
the structure and the instabilities of the disk. We will consider these effects and the influence to 
the spectral energy distribution of the radiation from the disk in the next paper. 


6. SUMMARY 

In this paper we have evaluated the density and temperature structures of a protoplanetary 
disk with the surface density of the minimum mass solar nebula, rotating with the Keplerian 
velocity around a central star, which has the typical parameters of T Tauri stars, and under 
the influences of the heating processes via viscosity and the star. The so-called a-model with 
a = 10“^ is used here for turbulent viscosity. Assuming the vertical hydrostatic equilibrium 
and the local radiative equilibrium, and utilizing the numerical calculation of the two-dimensional 
radiative transfer equation, we have obtained the following conclusions; 

1. The surface layer of the disk is heated and flared up via the irradiation from the inner hot 
disk rather than the central star in our disk model. 

2. At the outer region, x > lAU, the radiative heating from the inner disk dominates the 
viscous heating because the temperature at the surface region is high and the disk is geometrically 
thick, owing to the irradiation from the inner hot region of the disk. 

Making use of the resultant density and temperature profiles, we have examined the convective 
and the magnetorotational instabilities of the disk, which are expected to induce the turbulent mo¬ 
tion, and then, bring about the angular momentum and mass transfer, and/or the energy transport. 
As a result, we have found that; 

3. The whole disk is convectively stable since the dusty region is not heated by the viscous 
dissipation from the midplane but the radiative heating from the inner disk. 

4. Almost all of the disk is magnetorotationally unstable as the ionization degree is high 
enough that the stabilizations by ambipolar diffusion and ohmic dissipation are negligible, except 
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for the region near the midplane of 2AU < x < lOAU. 

Turbulent eddies, excited by magnetorotational instability, will trap dust grains to facilitate 
their growth under some conditions. We have discussed the size distribution of dust particles in 
the disk, which suggests that; 

5. The dust particles can grow to be cm-size in the surface layer region because of the high 
temperature distribution and the large disk height. These particles are expected to be exposed to 
the outside activities, such as x-ray radiations, the FU Ori outbursts, and so on, which may affect 
the dust grain evolutions. 

The influences of the evolution of turbulent motion, the dust size distribution, and the surface 
density distribution on the structure, the instabilities, and the spectral energy distribution of the 
disk should be studied in the next work. 
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clarity of our discussions. I would like to thank Dr. T. Nakamoto, Dr. T. Sano, and Dr. E.I. 
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and Dr. H. Kamaya for continuous encouragement. This work is financially supported by Research 
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Numerical Method 


A.l. Iterative Method to Calculate Temperature and Density Profiles 


At first we note how we evaluate iteratively the temperature and density distributions. In the 
following r*(r, 0) and / 0 *(x, z) represent the results of the f-th iterative calculation of the processes 
(II)-(V), while is the result of the j-th iteration of (a)-(f). The iterative methods are as 

follows; 

(I) Initially put the density and temperature profiles of optically thick accretion disk, z) 

and 0). They are derived in the conventional manner from the hydrostatic equation (1) in 

§2.1 and the energy transfer equation in the vertical direction via diffusion (e.g., Lin & Papaloizou 
1980, 1985), 

4c d{aT^) 


= (tT: 


3kkp dz 

where a is the radiative constant, defined with the Stefine-Boltzmann constant a and the light 
speed c as a = 4a/c. In order to solve equation (Al) we use a condition (9/4 )Sq!CsoD = crT^g, 
where Cgg = kT{x,z = 0)/m^ and = T{x,z = Zp). Zp is the photospheric disk height where 
the optical depth from the disk surface Zoo becomes 2/3, that is, K-^pdz = 2/3. Here 

kr is the Rosseland mean opacity described in §2.4. The condition means the energy balance in 
the vertical direction between the viscous heating at the midplane and the blackbody radiative 



- 13 - 


cooling at z = Zp (e.g., Meyer & Meyer-Hofmeister 1982). These initial density and temperature 
distributions are plotted and used in §3 and §4 as ‘non-irradiated disk’ (which means that only 
vertical energy transport is took into account), and compared with the ‘irradiated disk’ (obtained 
from our two-dimensional radiative transfer calculation below). 


(II) Evaluate the temperature distribution T*(^r,^0) by using z) and iteratively 

solving equations (7)-(14) as we mentioned in §2.3: (a) Calculate the intensity at each point 
(j)) by using tentative temperature distribution and integrating equa¬ 

tion (13) with the Short Characteristics method described in §A.2. (b) Substitute the intensity at 
a point Iy{r,Q-^ (f)) into equation (8) to obtain the mean intensity Jy{r,Q) at that point, (c) 
Evaluate the corrected temperature T*^(r,Q) at that point from equation of radiative equilibrium 
(7). (d) Estimate the temperature correction in context of the Newton-Raphson method (e.g., Press 
et al. 1985) as 


6P{r,e) 


{dAr/OTy 


{A^-^ - A^)/5T^-^ 


for A^(r, 0) > 0.1 and = 0 for A^(r, 0) < 0.1, where A:^(r, 0) = {r*-^ (r, 0) — 

T-^(r, 0)}/r-^ (r, 0). (e) Compute ST^^ry 0) for all grid points in the same manner, (f) Return to 
step (a) with the new tentative temperature r-^'''^(r, 0) = T-^ (r, 0) -|- ST^r, 0) at each point, and 
repeat these operations until |A^(^r,^0)| < 0.1 and |T'^(r, 0) — T'^“^(r, 0)|/T-^“^(r, 0) < 0.1 for 

(V0). 


(Ill) Making use of r*(^r,^0), compute the density distribution z) from equation (1) 

as we stated in §2.1. 


(IV) Estimate r*'’'^(r, 0) in the optically thick region by using p^i^x^ z) and solving equation 

(Al) as we stated in the process (I). This is because we can hardly get the adequate temperature 
from a wrong initial value of 0) in very optically thick region, owing to the approximate 

treatment of integral (13) and to the determination of the temperature correction ST^(r,Q) in 
(II). As we describe in §A.2, the Short Characteristics method to integrate equation (13) lets 
I,y{r,Q-, p,(f>) ByT{r,Q)], which leads to A:^(r, 0) ^ 0 and then STyr,Q) 0, for very 
optically thick limit. (This is also why we choose the initial density and temperature distributions 
as [I].) 

(V) Return to step (II) until \p'^{x,z)—p^~^{x,z)\/p'^~^{x,z) < 0.1 and |T*(r, 0)—T*“^(r, 0)|/T*“^(r, 0) < 
0.1 for (^xy z) and (^r,^0). 


A.2. Short Characteristics Method to Calculate Intensity 

Next we mention the numerical method to integrate equation (13) in order to get the specific 
intensity Since it costs vast CPU time to compute the integration from s' = 0 to s 

in all directions (p,, (p) and all frequencies ly, we here adopt an efficient and convenient algorithm, 
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the Short Characteristics method in the spherical coordinates (Dullemond & Turolla 2000; see also 
Mihalas, Auer, & Mihalas 1978, Olson & Kunasz 1987, Auer & Paletou 1994). Following them, 
we approximate to calculate the intensity at a point P under the assumption of neglecting the 
scattering (that is, only the black body radiation B^[T{r,Q)] is considered for source function at 
each point) as 


K, e '"''Iy{U-,ix,4))+UuBy{U)+pyBy{P) + duBy{D), (A2) 

where U and D represent the upstream and downstream points along a ray directed to {p, (p), and 
Ti, is the optical depth between the points U and P. The coefficients u^, Pu, and di, are functions of 
the optical depths from U to P and from P to D (see e.g., Dullemond &: Turolla 2000 for detailed 
determination of the coefficients). The intensity (A2) approaches Iy{P;p^(p) —> for 

optically thin limit, Ti, —> 0, and Iy{P;p^(p) —> By{P) for optically thick limit, Ty —> oo. Making 
use of this method, we can obtain 0;^ cp) systematically from the intensity at boundaries 

(see below) and the temperature at each point T(^r^ Q). 

In this calculation we use the spatial grid of (r^,©;), where k = 1,40 and I = 1,80. In 
order to resolve the inner region of the disk, the radial grid is put logarithmically as r^+i = pr^, 
where p is a constant, determined as ri = and r 4 o = rout are satished. We assign the inner 
radius as r\a = i?* = 2i?0 and the outer radius as rout = lOOAU. The azimuthal grid points are 
equally spaced in 0 < 0 < 7r/2. For the direction of the ray-path {p,<p), we put 32 points at the 
maximum for in 0 < ^ < 1 (see Dullemond &: Turolla 2000 in detail for choosing the ^-grid), 
and equally spaced 8 points for cp in 0 < (p < 7r/2. The frequency is divided logarithmically into 
32 points in 10^*^ < zz < 10^® Hz. The boundary conditions that we adopt here are as follows: 
Iy{r = ri, 0; ^ > 0, 0) = By{T = T* = 4000K) is taken for the intensity to the outward direction at 
the inner radius r = ri, and Iy{r = r 4 o, Q; p < 0,(p) = By {T = lOK) for the inward intensity at the 
outer radius r = r 4 o (the gaseous matter in molecular clouds is considered to have the equilibrium 
temperature of about T = lOK between heating by cosmic rays and cooling by molecular line 
emissions; e.g., Myers 1978). At 0 = 0i we take the symmetric condition against z-axis (0 = 0), 
and at 0 = 7r/2 we put Iy{r, Q = 7r/2] p, (p) = By {T = Tvis). We note that the boundary conditions 
at r = and 0 = 7r/2 reproduce the heating sources, the central star and the viscous dissipation 
at the midplane, respectively, that we described in §2.2 and §2.3. 
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Fig. 1.— The temperature {solid lines) and density {dotted lines) distributions in the x-z plane of 
(a) ‘irradiated disk’ (obtained from 2-D radiative transfer calculation) and (b) ‘non-irradiated disk’ 
(derived from only vertical energy transfer equation). The contour levels are T = 10,10^, 10^, 10^, 
and lO^K for the temperature, and p = 10“^®, 10“^®, 10“^^, 10“^^, and 10“®g cm“^ for the 
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Fig. 4.— The convectively {dotted stripe) and the magnetorotationally {solid stripe) unstable regions 
in the x-z plane of (a) ‘irradiated disk’ and (b) ‘non-irradiated disk’. The whole ‘irradiated disk’ 
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Fig. 5.— The maximum dust size distribution in the x-z plane of ‘irradiated disk’. The contour 
levels are Up^max = 10“^, 10“^, 1,10, and lO^cm. The thick line displays z = Zoo- The dust particles 
can grow to be cm-size in the surface layer, which are expected to be under the influences of the 
outside activities. 



